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Abstract. We review recent results on the nonlinear development of thermal insta- 
bility (TI) in the context of the turbulent atomic interstellar medium (ISM), in which 
correlated density and velocity fluctuations are present, as well as forces other than 
the thermal pressure gradient. First, we present a brief summary of the linear theory, 
remarking that, in the atomic ISM, the condensation mode is unstable but the wave 
mode is stable at small scales. Next, we revisit the growth of isolated entropy perturba- 
tions in initially unstable g function of the ratio of the cooling to the dynamical 
crossing times rj. The time for the dynamical transient state to subside ranges from 4 
to 30 Myr for initial density perturbations of 20% and sizes 3 to 75 pc. When <C 1, 
the condensation produces locally supersonic motions and a shock propagates off the 
condensation, bringing the surrounding medium out of thermal equilibrium. Third, we 
consider the evolution of velocity perturbations, maintained by a random forcing, repre- 
senting turbulent energy injection to the ISM from stellar sources. These perturbations 
correspond to the wave mode, and are stable at moderate amplitudes and small scales, 
as confirmed numerically. 

We then consider the behavior of magnetic pressure in turbulent regimes. Various 
observational and numerical results suggest that the magnetic pressure does not cor- 
relate well with density at low and intermediate densities. We propose that this is a 
consequence of the slow and fast modes of nonlinear MHD waves being characterized 
by different scalings of the magnetic field strength versus density. This lack of corre- 
lation suggests that, in fully turbulent regimes, the magnetic field may not be a very 
efficient source of pressure, and that polytropic descriptions of magnetic pressure are 
probably not adequate. 

Finally, we discuss simulations of the ISM (and resolution issues) tailored to in- 
vestigate the possible existence of significant amounts of gas in the "lukewarm" tem- 
perature range between the warm and cold stable phases. The mass fraction in this 
range increases, and the phase segregation decreases, as smaller scales are considered. 
We attribute this to two facts: the enhanced stability of moderate, adiabatic-like ve- 
locity fluctuations with ?7 3> 1 and the recycling of gas from the dense to the diffuse 
phase by stellar energy injection. Moreover, the magnetic field is not strongly turbulent 
there, possibly providing additional stability. We conclude by suggesting that the gas 
with unstable temperatures can be observationally distinguished through simultaneous 
determination of two of its thermodynamic variables. 
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1 Introduction 

The fact that the neutral atomic intersteUar medium (ISM) is most Ukely ther- 
mally bistable [ |68|p^ , |8^ has had a great impact on our picture of interstellar 
structure formation. Indeed, in two of the most influential models of the ISM 
to date, the famous two- and three-phase models of the ISM of Field, Gold- 
smith & Habing and McKee & Ostriker |Q, the concepts of thermal and 
pressure equilibria played a fundamental role, so that distinct phases (thermo- 
dynamic regimes with different density and temperature, but the same pressure) 
were predicted to coexist in pressure equilibrium. These phases correspond to 
stable thermal-equilibrium (i.e., heating-cooling balance) temperature regimes, 
and are separated by unstable regimes that, in those models, were therefore not 
expected to be present in the ISM. An opposite view was taken in the so-called 
time-dependent model of the ISM of Gerola, Kafatos & McCray |Q , which was 
presented as an alternative to the pressure equilibrium two-phase model, and 
which made radically different assumptions: a constant density in the presence 
of stochastic, local heating events that caused strong local fluctuations of pres- 
sure and temperature, because the cooling and recombination times are com- 
parable or shorter than the time between successive exposures of a given gas 
parcel to those heating events. This model predicted that significant amounts 
of gas should be in the unstable range, as they cooled after the transient heat- 
ing events. More recently, Lioure & Chieze [Q have considered models with a 
continuous recycling of gas among the various gas phases due to stellar energy 
injection, also concluding that significant amounts of gas should populate the 
unstable temperature range in the ISM. Note that the three-phase model 
did consider the existence of local fluctuations in the pressure, although it was 
still based on the premise of "rough pressure balance" . 

Nevertheless, both the equilibrium and the time-dependent models omitted 
a number of important aspects in the ISM budget. The multiphase equilibrium 
models essentially neglected the possibility of large pressure fiuctuations in the 
ISM. The time dependent model instead included this possibility as a funda- 
mental premise, but neglected the fact that such pressure fluctuations should 
induce motions, which should in general be turbulent (i.e., spanning a wide 
range of scales), and in turn cause strong density fluctuations [pOP]. Moreover, 
both the time-dependent and the three-phase models omitted other important 
agents of the ISM, such as magnetic fields, rotation, and cosmic rays. Elmegreen 
[ |l7||20|] performed a combined instability analysis including self-gravity, cooling 
and heating, and magnetic fields, but the effects of turbulence, which is an in- 
herently nonlinear phenomenon, can only be dealt with by means of numerical 
simulations of the gas dynamics in the Galactic disk in the presence of ther- 
mal instability (TI). The role of the turbulent motions may be crucial. In fact, 
realistic cloud/intercloud structure has been reported in models incorporating 
turbulence from stellar-like driving and cooling, but not necessarily a thermally 
bistable regime [p| jTl| , [l^ , |69t[78y66| , ^ . Similar results have been reported for pres- 
sureless (Burgers-like) models with stellar driving ||7^,|l^, and for simulations of 
interacting nonlinear MHD waves Thus, it is important to investigate the 
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role of TI in determining the distribution of the physical variables (density, tem- 
perature, velocity) of the flow, and, in particular, the degree to which phase 
segregation, as was proposed in the multiphase models, is realized, in the con- 
text of a turbulent ISM with multiple sources of turbulent energy at a variety 
of scales, such as stellar winds, supernova explosions, spiral arm passage, mag- 
netorotational instabilities etc., besides TI. 

Although the nonlinear development of TI has been studied extensively for 
decades now (e.g., |§|7|,|^,||,|l]j7§|l|,||,|§|,P,||,|§§), only recent work 
has started to investigate the interplay between TI and the turbulent nature 
of the ISM, such as, for example, the triggering of TI by external compressions 
( l3^ , |3^ , ^ ) and the possibility that the TI itself may contribute to the generation 
of turbulence in the ISM 

However, the fact that additional energy sources feeding ISM turbulence be- 
sides TI itself, such as stellar energy injection, or large-scale gravitational or 
magnetic instabilities, has additional implications. First, the very presence of 
strong motions implies that transport (advection) should be important, while 
traditionally conductive processes have received more attention (e.g., p3|j6||55t). 
Second, these transport processes may imply the existence of constant frac- 
tions of gas transiting through the unstable regime, and erase, to some extent, 
the phase segregation expected in multiphase models. These expectations are 
furthered by several observational studies (e.g., | |l^ , ^ , ^7|j25| , |3l|] ) that have sug- 
gested that the fraction of gas in the unstable range between the cold and warm 
phases of the atomic ISM is substantial. 

Another important issue to consider is the fact that the ISM is magne- 
tized, which suggests the possibility that turbulent magnetic pressure may sup- 
plement thermal pressure and somehow counteract TI. The condensation pro- 
cess in a magnetized medium has been studied by several workers (see, e.g. 
[ p3pl| , |6l| , p0p^ ; see also the references given in Q), concluding that, although 
condensation can be inhibited under some circumstances, it is in general possi- 
ble. However, those studies have not considered the case of TI developing in an 
externally-driven turbulent medium, except for [js^ . 

In this paper we review recent work and present new results concerning the 
interplay between TI and turbulence in the warm and cool ISM. First, we review 
the main aspects of the instability in Then, in §|[ we revisit the growth of 
isolated density fluctuations, focusing in particular on the late stages and the 
state of the gas surrounding the condensation. Next we discuss the development 
or suppression of growth in the presence of random velocity fluctuations, stress- 
ing that these probably constitute the most common way of inducing density 
fluctuations in the ISM. We then briefly discuss the nature of magnetic pressure 
on turbulent media (§||) and its dependence on density. In we discuss the 
role of TI in numerical models of the ISM in the presence of magnetic fields, 
the Coriolis force, modeled star formation, self-gravity and TI, aiming at deter- 
mining the fraction of unstable gas, and at interpreting the results in the light 
of the previous sections. We include an extensive discussion of numerical tests 
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to maximize the reliability of the results. Finally, we summarize the results and 
mention a number of implications in §|^. 

2 Review of the linear theory and the role of the ratio of 
cooling time to crossing time 

The ISM in general is subject to heating and cooling processes that cause it to 
have an effective thermal behavior very different from that of an ideal gas. Our 
current understanding of the physical processes responsible for the heating and 
cooling of the neutral atomic ISM can be found in Wolfire et al. . Their effect 
is customarily collected in a net cooling function per unit mass L — pA — F, 
where p is the gas mass density, A is the cooling rate, and F is the heating rate. 
In this paper we will use a piecewise power-law fit to the "standard" cooling 



curve of 1 82 , of the form 



yi = for T, <r<r,+i, (1) 

with the coefficients and exponents given in Table |l|, and a constant heating rate 
F = Fq, which is a reasonable approximation to the weak dependence F ~ 
across the density range of the atomic medium . The fit is obtained by first 
fitting a piecewise power law to the standard thermal equilibrium-pressure vs. 
density curve of [|2| to obtain the exponents of the cooling curve, and then 
determining the coefficients by equating the cooling rate to the constant heat- 
ing Fq = 0.015 erg s~^g^^. The values of the coefficients have units of erg 
s-ig-^cm^K-'^-'+i. 



Table 1. Cooling function parameters 



Interval 




rii/ cm 








15.0 


80.0 






(1,2) 






3.42 X 10^*^ 


2.13 




141. 


7.00 






(2,3) 






9.10 X 10^* 


1.00 




313. 


3.16 






(3,4) 






1.11 X 10^° 


0.565 




6101. 


0.59 






(4,5) 






2.00 X 10** 


3.67 




lO'^ 









" In units of erg s -"^g ^cm^K ''i.'+i. 
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Subject to these processes, and in the absence of magnetic fields, self-gravity 
and other agents, the dynamics of the gas are described by the equations 123,d5| 



dp 



pV • u = 0, 



dt 
du 

1 dP 



J — 1 dt 7^1 



PV ■ u + pC{p, T) - V • (ifVT) = 0, 



(2) 
(3) 

(4) 



where u is the fluid velocity, P is the thermal pressure, T is the temperature, 
7 is the heat capacity ratio, K — K{T) is the thermal conductivity, and the 
ideal-gas equation of state, P = pRT/ fi, is assumed, with R being the universal 
gas constant and p the mean molecular weight in units of the hydrogen mass. 
For clarity, we remind the reader that R = ks/mYi, where /cb is the Boltzmann 
constant, and tor is the mass of the Hydrogen atom. 

The linear analysis leading to TI was first performed in full in the classic 
paper of Field , and later generalized to the case of a flow in motion (sj,^ . 
Other useful, more recent presentations may be found in [ p"8||7^ , p5| | . The analysis, 
assuming perturbations proportional to exp(nt -|- ik • x) in all variables, yields 
the cubic dispersion relation 



Cv 



nc^k^ 



Np 
Cp 



ikK 



= 0, 



(5) 



where n is the growth rate, k is the wavenumber, c 
sound speed, and 



jkT/p is the adiabatic 



and 



R cpo 



(6) 



(7) 



pi-f- l)K 

As pointed out by Field, kx is the mean free path of the gas molecules (see also 
[^, p.35). 

The dispersion relation (||) has three roots, one of them being always real, and 
the other two being either a complex conjugate pair, or a pair of real numbers. 
There is instability whenever 5R(n) > 0, where denotes the real part. The 
positive real root corresponds to exponential growth without propagation, and 
is thus called the condensation mode. The pair of complex roots corresponds 
to oscillatory behavior, and is thus called the wave mode. This mode grows in 
amplitude (i.e., is overstahle) when the real part of those roots is positive. This 
nomenclature was extended by Field to the case when the roots are real, in which 
case he said the wave mode is overdamped (i.e., does not oscillate). 

The condensation mode is unstable if the so called "isobaric" criterion, namely 



Np < 0, 



(8) 
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is satisfied. However, when only this criterion is satisfied, the growth rate van- 
ishes as fc ^ 0. For the growth rate to remain finite at long wavelengths, it is 
necessary to also satisfy the "isochoric" criterion, which reads 

Np < 0. (9) 

The corresponding instability criterion for the wave mode is the so called 
"isentropic" (or "adiabatic") criterion, reading 

^-^>0, (10) 
Cp Cv 

where cp and cy are the specific heats at constant pressure and at constant 
volume, respectively. Note that cy = i?/ [(7 — For the coohng and heating 
functions adopted here, the isobaric, isochoric and isentropic criteria respectively 
imply /3 < 1, /? < and /3 < 1/(1 - 7) = -3/2. 

Since in the atomic ISM only the isobaric criterion is satisfied, and given 
the focus of this paper on that medium, we will concentrate on this case in 
what follows. It is worth recalling that isobaric and isochoric perturbations (ac- 
complished, for example, by setting up density fluctuations at constant pressure 
or temperature fluctuations at constant density, respectively) , are generically re- 
ferred to as enirop?/ preturbations, since they imply a variation of the ratio P/p"', 
which remains constant for reversible isentropic processes. We now discuss the 
evolution of entropy perturbations in detail. 



2.1 Entropy perturbations 

Figure shows the growth rate as a function of wavenumber for our fit to the 
"standard" cooling curve of |^ , eq. (^ and a realistic value of the conductivity 
K = Ko = IQ-^T^^^ erg cm^^ s'^ K'^ |l|, with To = 2400 K. In the case of 
the pure development of the instability, without any external forcing processes, 
there are three clearly distinct scale ranges for the wavelength A that arise from 
the presence of three characteristic time scales [ p5[ : the dynamical time, t^, in 
this case given by the sound crossing time, Tg — c/A; the cooling time, which for 
isobaric processes is given byR 



7/? 



ij-l)ti\Np\' 



(11) 



but which, to order of magnitude, is in general Tc ~ cyT / [pA)] and the con- 
ductive time, Tk — A^/kq, where — Kg/po is the thermal diffusivity. These 
time scales then define two characteristic length scales. First, the so called Field 
length Af = 27r/fcF, with fcp = (| A^p|/ko)^^^, at which r„ ^ Tc, and below 
which the growth of perturbations is inhibited by thermal conduction. Second, 

^ This expression is obtained by linearizing the energy equation and then comput- 
ing Tc = 5e/{dSe/dt) for an isobaric process. 
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the scale at which Td Tc, which wc shall denote Acq. In the remainder of the 
paper we shall use the notation rj = Tc/t^. The three scale ranges are then the 
"long" (A > Acq), "intermediate" (Acq > A > Af) and "short" (Af ^ A), where 
A = 27r/fc is the perturbation wavelength. For the long- wavelength range, 77 ^ 1, 
while for the intermediate- and short- wavelength cases, 3> 1. In what follows 
we will often discuss in terms of 77, as we consider that the relevant physical 
quantities involved are the cooling and sound-crossing time scales, even though 
it is customary in the literature to base the discussion on the perturbation wave- 
length. 

In the case of short- and intermediate-wavelength perturbations 0, the con- 
densation mode evolves nearly isobarically, as Tc ^ and thus sound waves 
have ample time to restore pressure equilibrium while the gas cools. This also 
has the consequence that, in the case of zero diffusivity, the growth rate asymp- 
totically approaches the cooling rate in the limit of large wavenumbers. In the 
presence of diffusivity, the growth rate decreases again in the short-wavelength 
range, due to the action of thermal conduction. These properties are illustrated 
in fig. |l|. 

In the opposite limit of very large wavelengths, the condensation mode ini- 
tially behaves isochorically (even though the isochoric criterion is not satisfied), 

^ For convenience, in the remainder of this paper, we will group intermediate- and 
small-scale perturbations into the small- wavelength (rj > 1) category. 




0.0 1 ... 1.. 

0-001 0-010 0,100 1.000 -0-000 100-000 

1/l[pc] 



Fig. 1. Growth rate as a function of wavenumber for our piecewise power-law fit to 
the "standard" cooling function of Wolfire et al. ^] and a constant heating rate. The 
wavenumber at which the growth rate becomes zero corresponds to the Field length. 
The points show the numerical growth rates in simulations of the li near growth of 
2.5% isobaric density perturbations in boxes of physical size 150 pc (see ^6.1) for several 
values of the mass diffusion coeflicient: ^ = {crosses); jj, — 0.001 (asterisks); ^ = 0.002 
(diamonds); /j, — 0.004 (triangles); /x = 0.008 (squares). For perturbations with / — 
1/0.05 = 18.75 pc (1/8 of the box size), which are in the transition regime between long 
and short wavelengths, growth rates are shown that have either zero (appropriate for 
the short-wavelength regime; points lying below the theoretical curve) or self-consistent 
velocity and pressure fluctuations (appropriate for the long-wavelength regime; points 
lying above the curve). 
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since 77 ^ 1, so that the coohng acts much more rapidly than the sound waves 
can travel accross the perturbation to restore pressure balance. This implies that 
large pressure gradients can be set up, which in turn can trigger strong motions 
that can become locally supersonic 54, |,^. Moreover, in this regime, the 
thermal pressure is given by the condition of thermal equilibrium, because the 
rapid cooling always allows its establishment. Figure ^ shows the equilibrium- 
pressure versus density for our piecewise cooling function (|l|) . The slope of this 
graph constitutes an effective polytropic exponent given by 7i,i+i = 1 — 
(c.f. eq. so that the pressure behaves as P oc pT^.^+i. In this figure, we de- 
note by pisob the density value within the cold phase that corresponds to the 
same pressure as that at the mean density. The instability under the isobaric 
mode is seen as the negative-slope range 0.6 cm~'^ ^ p ^ 3.2 cm~^ (equivalent 
in thermal equilibrium to the temperature range 300 ^ T ^ 6000 K). Note also 
the marginally stable, 7cff — behavior in the range 3.2 cm~'^ ^ p ^ 7.1 cm~^. 

In this same limit (long wavelengths), the growth rate asymptotically ap- 
proaches HI] 

/ \ 1/2 

niong = ±kc y~jf-J ' (12) 

which is of the order of the inverse of the sound crossing time, explaining its 
vanishing as fc (fig. |l|). Note that the situation is different when the isochoric 



4.0 




Fig. 2. Thermal-equilibrium piecewise polytropic behavior of the pressure for the cool- 
ing functions used in Papers I and IV. Pisob is the value of the density within the 
cold phase whose corresponding thermal-equilibrium pressure equals that of the mean 
density po- 
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criterion is satisfied, in which case, the growth rate remains roughly constant at 
over the long- and intermediate-wavelength ranges. 

2.2 Adiabatic perturbations 

The case of adiabatic perturbations is perhaps the most relevant for a turbulent 
ISM since the adiabatic condition implies that the density and pressure gradients 
have the same sign, a situation which is naturally accomplished by means of a 
compressive motion acting on time scales much shorter than the cooling time. In 
contrast, for entropy perturbations the density and pressure gradients can have 
opposite signs. We refer to the latter case as a "reversed" pressure gradient, 
while we say that a pressure gradient with the same sign as the density gradient 
is "regular". 

In the presence of velocity fluctuations, a new characteristic time scale ap- 
pears in the system, namely the bulk- velocity crossing time, = A/m, where u is 
the characteristic velocity of the perturbation. Thus, it is convenient to redefine 
the dynamical time scale as Td = min(rs, Tu), so that the definition rj — Tc/t^ can 
be preserved. Adiabatic perturbations with 3> 1 become unstable when the 
adiabatic instability criterion is satisfied [ p3[[7^ , and excite the wave mode 
of the instability, which in this case consists of nearly adiabatic sound waves 
with growth rate n = ike + {l/2)(Np/cp — Np/cy), where the imaginary part 
gives the propagation speed and the real part gives a modulation that grows 
only if the adiabatic criterion is satisfied (overstability) . However, at 77 < 1 (long 
wavelengths), the complex conjugate roots become real, and the wave mode be- 
comes a condensation mode, which is unstable whenever Np/{jNp) < For 
Np > 0, as is usually the case, this reduces to the isobaric criterion. Thus, in the 
atomic ISM, wave-like (or "velocity") perturbations are linearly unstable only for 
T] < 1, i.e, in the long-wavelength limit. For nonlinear perturbations, however, 
the density increase in the waves can accelerate the cooling and locally cause rj to 
become < 1, allowing the instability to proceed even in cases of perturbations of 
initial short-wavelength perturbations. A similar effect can occur even when the 
compression acts on the warm stable phase [ ^2p^ . We will discuss this further 
in§|. 

2.3 Entropy vs. adiabatic fluctuations 

It is now important to ask what kind of processes in the real world can generate 
the two kinds of fluctuations: entropy or adiabatic. In the absence of initial 
fluid motions, entropy fluctuations can only be produced by locally varying the 
cooling-to-heating ratio. On the other hand, if velocity fluctuations are used as 
the driver of the density fluctuations, as occurs in a turbulent medium, then the 
density fluctuations can either behave as entropy or as adiabatic fluctuations 
depending on scale and on the velocity amplitude. The production of entropy- 
like perturbations clearly requires that Tu 3> Tc, so that the thermal response 
of the flow proceeds under thermal equilibrium. At small scales, where 77 3> 1, 
this then implies that Tu 3> Tc ^ Ts, so that the motions have to be essentially 
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quasi-static. At large scales, where 77 ^ 1, we see that even supersonic motions 
can produce entropy-like fluctuations, as long as Tu,Ts ^ t^. The remaining 
possibility, i.e, Tu ~ Tg <C Tc, occurring for finite-amplitude velocity fluctuations 
at small scales, causes adiabatic-like perturbations. 

Another important distinction is that, for entropy perturbations, the motions 
are driven by the thermal pressure gradient generated by the instability, and 
tend to restore pressure equilibrium. Instead, in the case of externally-driven 
velocity perturbations, the motions drive the density and pressure gradients, 
which then feed back on the cooling and the motions themselves. Thus, the 
cause-effect relationship between the motions and the thermal pressure gradient 
are essentially reversed in the two cases. 



2.4 The magnetic case 

The linear instability analysis in the presence of a uniform magnetic field B was 
also studied by Field [p3| . Here we just briefly summarize his main results, and 
then discuss some recent work in the nonlinear regime. 

Qualitatively, Field concluded that the inclusion of the magnetic field should 
introduce three main modifications to the non-magnetic results. First, the wave 
mode splits into three modes, which correspond to the three modes of MHD 
waves: Alfven, slow and fast. Of these, the Alfven mode is "neutral", in the sense 
that it does not interact with the instability, since it is strictly non-compressive, 
while the fast and slow modes are governed by the same isentropic criterion. 
Second, the condensation mode is unaffected when the vector wavenumber k is 
parallel to B, while the field has a stabilizing effect when k is perpendicular to 
B. Finally, heat conduction is greatly reduced in the direction perpendicular to 
the field because of the spiraling of electrons between collisions. In summary, no 
major modifications to the overall picture were foreseen by Field, even though 
the dispersion relation changes from cubic to fifth-degree. These results were 
verified numerically by Goldsmith pO| . 



More recently, Loewenstein |50 has extended Field's linear analysis to the 
case of a stratified backgound medium, showing that condensation modes do ex- 
ist in cooling fiows, and that, over a certain wavenumber range, the presence of 
the magnetic field suppresses the damping of the instability due to conduction. 
In work more closely related to our focus in this paper, Hennebelle & Perault 
have investigated the role of a strong compression wave, interpreted as a 
turbulent velocity fluctuation, on triggering the nonlinear instability in the lin- 
early stable diffuse phase in a magnetized medium already segregated into two 
phases. These authors showed that the instability can indeed be triggered when 
the directions of the compression and of the initial magnetic field are oblique, 
giving the threshold values of the angle between them for condensation to oc- 
cur at various values of the compression Mach number and of the magnetic field 
strength. In the forthcoming sections we will first describe results concerning the 
development of TI from an initially unstable medium under both quiescent and 
turbulent conditions, in the non-magnetic case. Next we discuss the nature of 
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magnetic pressure in turbulent media, and finally we will consider more complete 
models of the ISM incorporating all of these agents and processes. 

3 Nonlinear evolution of entropy perturbations 

As a first step in our discussion of dynamical aspects of the development of TI, 
in this section we discuss the conditions necessary for the development of large 
velocities (possibly supersonic) and shocks during the spontaneous (i.e., in the 
absence of external triggers) condensation process of entropy perturbations in 
the unstable atomic ISM, as a function of the parameter rj. A detailed review 
of the nature of the shocks has been presented in js^ . We are also interested in 
the duration of the dynamic phase. Studies dealing with the development of su- 
personic motions and shocks as a consequence of the condensation process from 
the unstable regime have mostly focused on the regimes of proto-galaxy-cluster, 
proto-galactic, and proto-globular cluster clouds (e.g., p^ , [70| , p8|j59| ,|8|, ^ , ^ ) . In 
particular, Sasorov [f70| pointed out that the development of TI in three dimen- 
sions should give rise to flattened structures, similar to the "pancakes" formed 
by gravitational contraction in the cosmological flow. 

In the context of the ISM, the nonlinear development of isolated entropy 
fluctuations was initially studied by Goldsmith |Q and Schwartz, McCray & 
Stein Q, who found that the condensation of small-scale (77 > 1) perturbations 
occurred on time scales of ~ 1-10 Myr, and produced clouds of densities 100 x 
larger than their initial values. Goldsmith also considered the case of large-scale 
perturbations, finding the development of transonic velocities, and that the time 
required for reaching a true steady state is much longer than the time of "initial 
collapse" , not being reached by any of his simulations. However, these works 
were performed at very low resolutions, and did not discuss the state of the 
surrounding gas in much detail. 

More recently, Burkert & Lin |^ have considered a cooling-only medium (i.e., 
without background heating), and suggested that a special clump scale can be 
selected by the following mechanism. In the case of a globally cooling medium 
which eventually exits the thermally unstable range in roughly one cooling time, 
large-scale (ry < 1) fluctuations that cool isochorically do not change their den- 
sity appreciably before exiting the unstable range, so that, after they do, the 
pressure gradient becomes regular again, and the perturbation is erased. On the 
other hand, small-scale fluctuations can eventually reach the regime of isochoric 
cooling with rj < 1 as their density and local cooling rate increase, and then 
stop growing, except if they reach this transition stage already with nonlinear 
amplitudes, in which case advection overtakes the pressure gradient in promot- 
ing the compression, which then proceeds at an accelerated rate. Thus, Burkert 
& Lin suggested that the latter fluctuations are the ones that dominate the 
fragmentation of a large cloud into clumps, determining the clump properties. 

In the presence of both cooling and heating, Sanchez-Salcedo et al. |Q (here- 
after Paper I) have recently investigated the evolution of perturbations as a 
function of size (or, equivalently, 77), focusing on the magnitudes of the veloci- 
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ties that develop, the time scales for reaching a relatively quiescent stage, and 
the state of the gas surrounding the condensation after the latter has reached 
a quasi-stationary state. To this end, Paper I performed one-dimensional (ID), 
high-resolution numerical simulations of the evolution of Gaussian-shaped per- 
turbations in the absence of any other physical processes. The simulations use 
up to 7000 grid points, and non-uniform grid spacing, with central grid point 
density ^ 100 x that at the edges, in order to maximize the resolution at the 
condensation center. The simulations solve the gas dynamic equations in the 
presence of heating and cooling parameterized as in eq. (jl|) , reading 

-Plnp _ ^ ..ON 

Dt " dx^ p^'dx^' ^ ' 

+ / + (14) 
Dt p ox p ox 

Ds 

T— ^ r- pA + 2vS'^ + diffusion term, (15) 

where D/Dt — d/dt + ud/dx is the Lagrangian derivative, S — {2/3)du/dx is 
the generalized strain tensor in the ID case, s is the entropy per unit mass, / 
is the random forcing (described and used in and the other quantities have 
their usual meanings. A shock-capturing viscosity of the form 

= 1^0 + Cy8x^ max(0, —V • u) (16) 

is used, with Cy a constant. The last term in eq. (|l^) is an artificial mass diffu- 
sion term, necessary for smoothing excessively large density gradients. The term 
2z/5'^ in eq. ([l^) is the viscous heating of the gas. The diffusion term in the 
same equation has the same form as the mass diffusion term, and is included to 
guarantee that any mass redistribution occurs with the corresponding entropy 
redistribution. The coefficients /x of these terms are maintained at very low val- 
ues so that these diffusivities are comparable to the numerical diffusivity, giving 
a diffusive scale of size 2-3 grid zones. All the results reported in this section 
have been subjected to convergence tests to guarantee that they are not altered 
by changing the resolution (see Paper I). 

The simulations start at a density of 1 cm^'^, roughly the mean ISM density 



in our galaxy, as pointed out in § |2.4| from , which lies in the unstable range 
of the cooling curve. The equilibrium temperature at this density is T w 2400 
K. The scale Aeq at which the cooling and sound crossing times are equal is 
Acq ~ 10 PC, while the Field length under these conditions is Ap ^ 0.015 pc (cf. 
Qand fig.|l|). 



We first consider the case of small-scale {r\ > 1) entropy perturbations, as 
they constitute the paradigm of cloudlet (i.e., small clouds of sizes ^ 1 pc and 
densities - 50 cm'^) formation by TI in the ISM (e.g., ||,|o|j74|j6|,||l ) . To this 



end, we have performed a simulation, labeled DENS, of the evolution of a Gaus- 
sian density perturbation of 20% amplitude and a full width at half maximum 
(FWHM) of 3 pc. Figure ^ shows the density, pressure and velocity profiles of 



Thermal Instability and Magnetic Pressure in the Turbulent ISM 



13 



the cloud at various times until the time when a "cloud" has formed and the 
accretion process has mostly subsided. Note the logarithmic a:-axis, where x is 
the distance to the center of the cloud. It is seen that the evolution is indeed 
quasi-isobaric, with variations in the pressure of less than 2%, and local Mach 
numbers which do not exceed 0.2. By t ~ 4.2 Myr, the condensation has es- 
sentially completed its evolution, and reached the pressure-equilibrium density, 
Pisob- Figure |^ shows the evolution of this run on the P-p diagram. The quasi- 
isobaric nature of the condensation is clearly seen, especially at the final time, 
at which all points with densities higher than the mean have almost exactly the 
same pressure. 

Note, however, that in figs. || and ^ a population of points is still seen to 
continue flowing onto the condensation, and as it does, it necessarily remains 
in the "unstable" range. In fact, the mass of this gas amounts to ^ 8 times 
the mass in the condensation (counting only the gas reached by the rarefaction 
wave). In fig. |[ it is seen that this occurs in the region <L log \x/ (pc)| ^ — 1.5. 
The accretion and evacuation of the unstable gas will take very long times to 
complete, because the reservoir of unstable gas outside the cloud is very large 

90% of the total mass), in agreement with the remark by Goldsmith jSOt . 
Moreover, note that the inflowing gas is not truly unstable, as it does not lie on 
the equilibrium curve anymore. Instead, at t = 4.2 Myr, the density and pressure 
gradients have the same sign throughout this region. Thus this gas does not have 
a tendency to fragment any further, even though its density is in the "unstable" 



(°) 



1 1"" 


'"1 1 

t=4.2 

f 




t=3.4 




t=2.1 






1 i..<. 


t=0 


-1 
log Ix 


-2 -3 - 




E 0.4 



log Ixl [pc] 




Fig. 3. Time development of run DEN3. The density, pressure and velocity profiles 
at different times are plotted in panels (a), (b) and (c), respectively. The times corre- 
sponding to each line type are indicated in frame (a) in Myr, and the same line labeling 
is used in frames (b) and (c). 
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range. It can be said that this gas is flowing because of TI, but once it is doing 
so it has no tendency to fragment any further. 

Finally, note that the cloud formation time is not very short, and is signif- 
icantly sensitive to the initial amplitude. Simulations with an initial amplitude 
of 10% require ~ 5.5 Myr to complete the condensation. This time is compara- 
ble to the mean time between successive exposures to passing shock fronts from 
supernova remnants and superbubbles |^ , so that the condensations may have 
their growth interrupted by external perturbations, as is the case in 

Let us now consider the opposite case of a large-scale entropy perturbation, 
with FWHM=75 pc and ry ^ 0.04, in a box of 250 pc. We refer to this simulation 
as run DEN75. Its evolution is shown in figs. || (density, pressure and velocity pro- 
files) and^ (evolution on the P-p plane). In this case, the evolution is significantly 
different. As dictated by the smaller growth rates of larger-scale perturbations, 
run DEN75 requires 30 Myr to complete the formation of a cloud, but moreover, 
throughout the first part of its development, the condensation proceeds along 
the thermal equilibrium curve (see the first four panels of fig. developing 
locally supersonic velocities in the process (maximum Mach number ^ 1.2) that 
cause a strong overshoot. Thus, this condensation transiently reaches densities 
~ 55/9isob- At the time of maximum compression, a strong shock is produced 
at the cloud boundary that propagates outwards from it. This shock weakens 
quickly as it propagates into the low density medium, but it has the important 
effect of heating the still-infalling gas, bringing it out of thermal equilibrium 
and closer to isobaric conditions (cf. last two panels of fig. ||), and reversing the 
velocity gradient. This shock is located at the peak of the velocity for the various 
times shown in panel (/) of fig. |[ and is seen also in the pressure (see panel (e) 
of the same figure) . By the end of the simulation, the accretion ram pressure is 
still high enough that the condensation has p ~ Spisob, and this value decreases 




Fig. 4. Time development of run DENS in P-p phase space, at the times indicated 
in each frame. At t = 4.2 Myr, although the cloud has already formed, a substantial 
fraction of the points in the simulation are still traversing the unstable range, albeit in 
a nearly isobaric regime. 




Fig. 5. Same as fig. ^ but for run DEN75. Tlie first part of tfie evolution is sliown 
in tfie upper frames, and tfie rest in tfie lower frames. Note the formation of a shock 
shortly after t = 22.6 (frame c), which then propagates outwards from the cloud. At 
the same time, the density overshoots to over 55pisob (frame d). After the formation 
of the cloud, the density relaxes to a value ~ 2.5pisob, due to the ram pressure of the 
still infalling gas. 



extremely slowly with time. Again, as in run DENS, the infaUing gas is mostly 
in the "unstable" density range, in this case with a mass of almost twice that in 
the cloud. Within the infalling region, the pressure and density gradients have 
the same sign, so this gas again has no further tendency to fragment. We have 
found from other simulations that the qualitative behavior of run DEN75 occurs 
down to initial fluctuations with FWHM=15 pc. 

From the evolution of these two simulations, we conclude that large-scale 
(77 1) entropy perturbations have such a dynamic evolution that their fi- 
nal central density and pressure are larger than those corresponding to plain 
thermal-pressure equilibrium with the diffuse phase, and moreover require such 
long times to evolve (over 20 Myr to the occurrence of the large density over- 



16 Vazquez- Semadeni et al. 




Fig. 6. Time (icvclopmcnt of run DEN75 in P p phase space. Before the shock for- 
mation, the evolution proceeds along the thermal equilibrium curve. Subsequently, the 
outwards-propagating shock brings the outside medium out of thermal equilibrium, 
and restores nearly pressure balance. Thus, the infalling gas is traversing the density 
"unstable" range, but in nearly isobaric (inertial) conditions, rather than along the 
thermal equilibrium curve. 



shoot) , that they arc unlikely to complete their evolution before being disrupted 
by other perturbations in the real ISM, such as passing shock waves, or simply, 
general turbulent fluctuations. Small-scale entropy perturbations (77 ^ 1), on the 
other hand, adhere better to the paradigm of forming near pressure-equilibrium 
condensations, although we have seen that a significant fraction of the mass still 
lies in the unstable range after the cloud has formed, and is accreting onto the 
condensation, causing the presence of (weak) accretion fronts (rather than con- 
tact discontinuities) at the cloud boundaries that only subside asymptotically 
in time. Since the evacuation of the low density regions must proceed in times 
of order of the sound crossing time, the final fraction of mass in the "unstable" 
density range, in the more realistic case of multiple fluctuations, should depend 
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on their number. In fact, we have performed simulations with a full spectrum of 
initial fluctuations, and in those cases the final unstable fraction may be much 
lower, although still times ^ 15 (respectively, 8) Myr are required to evacuate 
the unstable range when the minimum perturbation size is 12.5 (respectively, 
1.25) pc. More importantly, however, small-scale perturbations behave very dif- 
ferently when they are quasi-adiabatic rather than quasi- isobaric, as we discuss 
in the next section. 



4 The case of velocity fluctuations 



As mentioned in § p.2| , velocity fluctuations are likely to be the most represen- 
tative of the actual situation in the turbulent ISM, because in a continuum any 
density fluctuation must originate from compressive or expansive motions. Such 
motions are readily available in a compressibly turbulent medium. When the 
cooling time is long {r] » 1), these compressions/rarefactions heat/cool the gas 
adiabatically, and the perturbations are then isentropic, which, as discussed in 



§2.2, are stable in the atomic ISM in the short-wavelength limit. As also men- 
tioned in that section, in the case of velocity fluctuations, the dynamical time in 
rj = Tc/ra is given by Td = minjrs, Tu}, where is the turbulent crossing time, 
which is in general also a scale-dependent quantity. 

Several studies |^,^,|3| have investigated the possibility of triggering TI 
through the nonlinear compression, either by strong shocks or large-scale, large- 
amplitude compressive waves of the warm stable phase, concluding in general 
that triggering TI off the stable phase is possible for strong enough compressions, 
with the possibility of even forming molecular hydrogen in the collapsed region 
However, these studies have assumed that the gas has already previously 
segregated into phases. It is our interest now to discuss the extent to which 
such segregation can be achieved, starting from unstable conditions. Therefore, 
in this section we describe the evolution of a thermally unstable medium (with 
respect to the isobaric criterion) subject to random velocity forcing, as originally 
presented in Paper I. 

We take uniform-density initial conditions, and apply a random forcing / of 
the form 

/(x, t) = Re [N exp (ifc (t) x + ic/) (t))] , (17) 

where k{t) is a time dependent wavenumber, (j}{t) is the phase and x the position. 
Following 0, we take N = /qCs {k{t)cs/dt)^^^ , where /o is a constant factor and 
6t is the length of the timestep. The values oik/ {2tt/1), where I is the box length, 
and of (j) are selected at each timestep randomly in the ranges [3, 10] and [0, 27r], 
respectively. The positive exponent (1/2) in N implies that strongest forcing 
occurs at the highest wavenumbers of the forced range, so that the energy- 
injection scale scale is Ai = I/IO. We have chosen this forcing for two main 
reasons. One, it mimics the small-scale stellar forcing acting in the ISM, and, 
two, it allows us to maintain the desired rms Mach numbers at the (small) scales 
of interest, since this is difficult to achieve with pure large-scale forcing. We do 
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not consider decaying-turbulcncc situations, as we are ultimately interested in 
models of the ISM, which is subject to continued energy-injection processes. 

We have performed nearly 20 simulations varying the box size and the scale 
and amplitude of the forcing (see Paper I for details). From them we conclude 
that, for the average conditions of the ISM, the presence of turbulent motions 
with small enough sizes (~ 0.3 pc) and moderate amplitudes (A^rms ~ 0.3) such 
that 77 is maintained above unity, condensations do not appear. We understand 
this as a consequence of the turbulent crossing time becoming shorter than the 
growth time of the condensations, allowing the turbulent fluctuations to both dis- 
rupt the incipient condensations and to more than compensate cooling through 
the heating from shocks and adiabatic compression; i.e., the perturbations be- 
come adiabatic, and therefore stable according to the linear analysis in the limit 
77 > 1. 

In general, we conjecture that the presence of velocity fluctuations in the 
ISM, even if unable to completely suppress the development of condensations, 
may increase the fraction of gas in the "unstable" temperature range. Note also 
that for small-scale velocity fluctuations, the evolution is not along the thermal- 
equilibrium curve, but rather intermediate between adiabatic and isobaric. Thus, 
the density range determined by the "unstable" temperature range under these 
conditions does not exactly coincide with the unstable density range under ther- 
mal equilibrium, as given in §^ and fig. |^. An interesting possible application for 
observationally measuring the actual thermodynamic state of the atomic ISM is 
mentioned in SfH. 



5 The magnetic pressure in turbulent media 

We now make a pause in the discussion of the thermal instability and consider 
the character of the magnetic pressure in turbulent flows, in order to assess the 
possibility that it may supplement the thermal pressure in the ISM, and thus 
contribute to weaken the effects of TI. Note that in this section we make no at- 
tempt to discuss TI in the presence of a magnetic field. This has been discussed 
by a number of authors (e.g., [^3|j6]] , |50| , ^ ). Instead, here we investigate the na- 
ture of magnetic pressure in fully turbulent compressible, magnetized isothermal 
flows. Several works have considered this regime as well, both numerically (see, 
e.g., the reviews by Mac Low, Ostriker and Nordlund in this volume, and refer- 
ences therein) , and theoretically |4£]. In particular, the numerical simulations of 
references and (see also [36| for the nonisothermal case) reported a lack 
of correlation between the density and the magnetic pressure, B^, where B is the 
magnetic field strength, at low and intermediate densities in cases in which the 
magnetic /3 parameter^, equal to the ratio of thermal to magnetic pressure, is 
^ 1. Moreover, recent observational data (see, e.g., Crutcher, Heiles & Troland, 
this volume) suggest a similar lack of correlation at densities below ~ 1000 cm~^ 
in the ISM. Paper II has attempted to understand the origin of this decorrelation 



3 



The magnetic /3 should not be confused with the exponent Pij of eq. ^ 
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in terms of the so-called "simple" MHD waves, and discussed its implications on 
the role of as a pressure. In this section we briefly summarize the results of 
that paper. 

We consider isothermal MHD flows in "1+2/3" dimensions, or slab geometry. 
The direction x is referred to as the direction of wave propagation. In this setup, 
bx, the field component along x, is constant. We denote by b the magnitude of the 
field component perpendicular to x. The initial, uniform magnetic field is chosen 
to lie in the {x, z) plane, at an angle 6 from the x axis, so that bx = cos 9 at all 
times. The treatment in this section is entirely in non-dimensional units, so that 
the parameters characterizing the flow are the sonic and Alfvenic Mach numbers 
of the velocity unit, denoted Mg and Ma, respectively, and the propagation angle 

0. The plasma beta is then P = M\/Ml^_^ 

"Simple" MHD waves (see, e.g., ||47| , |5^ ) are finite-amplitude solutions of the 
equations, characterized by the property that all variables can be expressed as 
wave "profiles", i.e., as a function of a single one of them (say, the density) as in 
the case of linear MHD waves. The same well-known modes of the linear case, 

1. e., Alfven, fast, and slow, exist in the case of simple waves. Only the latter two 
are associated with the density fluctuations. The propagation velocities of the 
modes are given by p7 52 



2MIp 



{B^ + /3p) 1 ± 



Afiblp 



(18) 
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Fig. 7. Magnetic pressure-density correlation, indicated by the two-dimensional his- 
togram of points in log-log coordinates, for a simulation with a magnetic field perpen- 
dicular to the direction of propagation (i.e., cos 6 = 0), and forcing parallel to this 
direction. This configuration allows only the existence of the fast mode of nonlinear 
MHD waves. The run has an rms field fluctuation SB/ B = 0.62 and 5p/ p ~ 0.62. The 
rms Alfvenic Mach number Ma = 5.2. The magnetic pressure is seen to scale as p^. 
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and 

where v± denotes the speed of the fast (+) and slow (— ) modes, and va, that of 
the Alfvcn mode. — + is the total field strength. 

After manipulating the equations to obtain the wave profiles, one finds, in 
particular, for the dependence of the field with density, 

d b^ d , o o ^ V 

In the limit when 4/36^p ^ {B^ + (3p)'^, this equation can be simplified and 
integrated using ([l8| ) to give the dependences for the fast and slow modes as 

b^ - 2(3p (slow mode) (21) 
B'^ oc p2 (fast mode), (22) 

where C is a constant. These equations essentially give the behavior of magnetic 
pressure with density for the two modes in the limit mentioned above. This 
condition is generally satisfied, except when j3p « 6^ and simultaneously b^ <^ b^, 
i.e., when b^ is not too small, for Pp of order unity and small field distortions. 

Several points are noteworthy about (|2^) and (|2|): 1) The /3 weighting in the 
second term of the RHS of (|2^) implies that at low /3, the slow mode produces 
large density fluctuations even when the field fiuctuations are small. 2) The 
total pressure, Ptot = + is rougly constant in the slow mode. 3) Most 
importantly, the pressures from the two modes depend very differently on density. 
One can then expect that, in the large fluctuation amplitude case (i.e., the fully 
nonlinear regime), the particular value of the magnetic pressure of a fluid parcel 
will not be uniquely determined by its density, but instead, that it will depend on 
the detailed history of how the density fluctuation was arrived at, causing a lack 
of correlation between the magnetic pressure and the density. 

The latter suggestions have been tested in Paper II by means of numerical 
simulations with random forcing (actually, an acceleration) applied on wavenum- 
bers 1-19 to all three velocity components or to only the perpendicular ones. 
Choosing the direction of the forcing and of the uniform magnetic field allows 
us to highlight either one of the slow and fast modes. Figure shows the b^- 
p correlation by means of isocontours of the two-dimensional histogram in the 
log(6^/2M^)-log(p) plane of the points in a simulation with 4096 grid points, 
forcing applied on all three velocity components, and the initial magnetic field 
perpendicular to the direction of propagation. This is a case in which only the 
fast mode exists, consisting of a pure magnetosonic wave, and the correlation 
exhibits the well known p^ behavior of magnetic pressure in this case. This result 
holds independently of the value of . 

Figure ||, on the other hand, shows the correlation for a run at the same 
resolution with the forcing perpendicular to the propagation direction, and the 
magnetic field almost parallel to the forcing {b^ — cos9 — 0.1), at low Alfvenic 
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Mach number (Ma ~ 0T5). In this case, the density fluctuation production is 
dominated by the slow mode because fi ^ 0.007 (c.f. (^l|)). The near constancy 
of the total pressure in this case is evident in this figure. In this simulation, 
one observes large oscillating density clumps which do not merge nor pass each 
other, mostly anti-correlated with b^. However, for this same field configuration, 
as Ma is increased, the field becomes more strongly distorted, and the fast mode 
starts acting on the perturbed field. The fluctuations become more random and 
superpose each other, with the presence of fast shocks and a correlation between 
p and b^. The result is that at high Ma both modes are actively producing 
density fluctuations, and the correlation between magnetic pressure and the 
density is lost, as seen in fig. ^. The situation can be idealized by assuming that 
the density of each fluid parcel is arrived at through a random sequence of slow 
and fast waves, which is different for each parcel. 

The case of parallel propagation is also interesting to mention as it illus- 
trates the complexity of this problem. For example, the picture of non-interacting 
clumps observed at small Ma (i.e. in presence of weak field distortions) and for 
large angles is also observed for parallel propagation in the case of large Ma- 
But in that case the clumps form inside strong slow shocks. The magnetic field 
intensity is the weakest inside the clumps which cannot merge due to the high 
external magnetic pressure. In the general 3D case, we expect all angles between 
the magnetic field and the propagation direction to be present, and therefore a 
representative case would be one in which the field is at 45° from the propaga- 
tion direction. In this case, we recover the trend of a roughly constant magnetic 
pressure at small Ma and an increased scatter between the magnetic pressure 
and density as Ma is increased (not shown). It is found that the level of den- 
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Fig. 8. Magnetic pressure-density correlation for a simulation with cos 6* — 0.1, and 
forcing perpendicular to the propagation direction, with 5B/B = 0.32 and Ma = 0.48. 
This configuration highlights the slow mode. The magnetic pressure is seen to remain 
virtually constant. 
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sity fluctuations greatly depends on the dominant mode, has no specific relation 

with (3 and is actually the largest when B is strong (and thus slightly perturbed). 
Density PDF is close to a log-normal when and p are not correlated, and 
shows an excess at small density when slow waves dominate. 

We conclude from this section that the scatter between magnetic pressure 
and density found in simulations and in observational data can be understood 
in terms of the different dependence of the magnetic intensity on density for 
the slow and fast modes of simple nonlinear MHD waves, and of the random 
sequence of these that a fluid parcel experiences as it evolves in a fully turbulent 
regime. Moreover, a number of important implications emerge from this lack of 
correlation. First, it suggests that modeling magnetic pressure by means of a 
polytropic dependence on density may not be adequate in the fully turbulent 
case. This relation is seen to apply when scale separation is preserved between 
small-scale Alfvn waves and large-scale density perturbations. Second, the mag- 
netic "pressure" does not really act as a pressure, as it does not behave as a 
restoring force in general. Instead, it acts more as a random forcing. Third, the 
latter point suggests that magnetic pressure should not be effective as a sub- 
stitute for thermal pressure when the latter behaves "in reverse" in thermally 
unstable situations. A flnal note is that the isothermality of the flows considered 
here is of no relevance to the results, which are thus expected to apply equally 
to non-isothermal flows. 




Fig. 9. Magnetic pressure-density correlation for a simulation with cos 9 = 0.1, and 
forcing in the three directions. This configuration also highlights the slow mode, but 

increases the field distortions {5B/B — 2.17 and Ma ~ 7.29). The magnetic pressure 
is seen to decorrelate from the density, although it appears to be bounded above and 
below by the slow and fast mode dependences, respectively. 
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6 TI in models of the ISM 

In the previous sections we have discussed the nonlinear development of fluctua- 
tions of various kinds and sizes in the presence of TI under the isobaric criterion, 
and the nature of magnetic pressure in turbulent media. With this background 
we can now proceed to discuss the behavior of numerical models of the ISM 
incorporating the magnetic field, self-gravity, rotation, shear, and stellar-like 
(localized) energy injection, in the presence of isobaric TI. 

As mentioned in the introduction, the classic two- and three-phase models of 
the ISM were based on the principles of thermal and pressure equilibrium, and 
thus did not predict the existence of significant amounts of gas in-between the 
phases. Instead, the time-dependent model [2^ did. Furthermore, observations 
do not clearly support a sharp phase segregation either; instead, they often 
have found evidence of significant amounts of "lukewarm"_gas at temperatures 
intermediate between those of the cold and warm phases (l|||,||,|5|j3|. Thus, 
it is important to perform numerical simulations that quantify this fraction and 
that allow us to determine whether sharp phase segregation is expected in the 
ISM, or whether it is more likely a continuum. 

Numerical simulations of the ISM including radiative heating and cooling 
as well as stellar-like energy injection have been performed by a number of au- 
thors, starting from the pioneering work of Bania & Lyon and continuing 
with references [0jl|,||,^j66H2|,^,|6|j4|,|l|,|^j27|,|^,||. These works have in- 



cluded different amounts of physics in the simulations, such as the magnetic 
field, self-gravity, galactic disk rotation, etc. However, the role of TI had been 
discussed only in passing until recently. Bania & Lyon, using two-dimensional 
(2D) simulations of a 180-pc square region on the Galactic plane, at a resolu- 
tion of 40 X 40 pixels, including randomly-positioned stellar-like energy sources, 
pointed out that the inclusion or not of a thermally unstable range had little 
effect on the resulting structure. Elmegreen |^ pointed out that high-resolution 
ID MHD simulations forced with nonlinear magnetic waves can form hierarchi- 
cal cloud/intercloud structure both if the cooling functions used have a single or 
double stable equilibria. The 2D non-magnetic simulations of Gerritsen & Icke 
[ p9| , and the 3D MHD ones of Korpi et al. |Q , although not specifically aimed 
at this issue, already pointed towards the existence of unstable gas. Vazquez- 
Semadeni, Gazol & Scalo |^ (hereafter Paper III) first discussed the interaction 
of TI with the turbulent motions in the ISM, suggesting that at large scales (sim- 
ulation size of 1 kpc) the signature of TI in the mass density histogram is erased 
when small-scale forcing mimicking expanding HII regions placed at the den- 
sity peaks is used in low-resolution (128 x 128) 2D simulations including the 
magnetic field, self-gravity and rotation. Gazol et al. (hereafter Paper IV) 
then reported that the temperature histogram in similar simulations at higher 
resolution (512 x 512) contains roughly half the mass at unstable temperatures. 
Recently, Mac Low et al. ^l[ (see also Mac Low, this volume) have presented 
3D MHD simulations spanning the whole range of temperatures existing in the 
ISM, including the hot (T Si 10^ K) gas, showing that the effect of supernovae 
is to introduce large fluctuations (by 2-3 orders of magnitude) in the thermal 
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pressure of the ISM, in contrast with the multiphase models of the ISM p^JS^ , 
and more in agreement with the time-dependent ones [ p8[ , and with the recent 
observational study of Jenkins & Tripp p6[ . Finally, in an analytical treatment 
of the effective equation of state of the ISM and the power spectrum of the en- 
ergy sources, Norman & Ferrara |Q concluded that the traditional multiphase 
description should be replaced by a "continuum of phases" . 

A related line of study has been that taken by Koyama & Inutsuka [Q 
and Kritsuk & Norman |46| . Both of these groups have recently considered the 
generation of turbulence in flows in more than one dimension, as a consequence 
of the nonlinear development of TI. Koyama & Inutsuka considered a shock- 
compressed layer in 2D between a fast flow with diffuse-gas properties {n = 0.6 
cm~^, T = 6000 K) and a hot gas region, showing that the layer fragments 
into small cloudlets that have supersonic velocity dispersions with respect to the 
warm medium in which they are embedded, and coalesce to form larger units. 
Kritsuk & Norman considered the 3D development of the instability alone, very 
far from thermal equilibrium (n = 1 cm^'^, T = 2 x 10^ K), so that the gas 
is initially unstable under the isochoric mode and the cooling times are very 
short (~ 0.3 Myr). They again found that the development of TI generates 
turbulence, in which roughly 15% of the mass is in the unstable regime. They 
point out, however, that this turbulence is decaying, because it has no other 
energy sources than the development of TI itself. 

However, none of the works reporting a fraction of unstable gas have discussed 
in detail the possible suppression of the instability by numerical limitations, 
raising a concern that perhaps the presence of unstable gas, and therefore the 
lack of sharp phase transitions in the simulations, are numerical artifacts. In the 
remainder of this section we first present new results concerning the numerical 
issues in detail, in order to assess the validity of the ISM simulations and interpret 
their results. We then present new simulations that, based on the numerical 
considerations, provide reasonable evidence that significant amounts of gas at 
unstable temperatures should be expected in the atomic ISM. 

The simulations presented in this section have been performed using a pseudo- 
spectral scheme, described in detail in [^, to solve the MHD equations in the 
presence of heating, cooling, stellar-like energy sources, and self-gravity, using a 
hyperviscosity scheme and including a mass diffusion term with coefficient /i (cf. 
eq. |l^) and a thermal diffusion term with a constant coefficient K = Kq (cf. eq. 

In non-dimensional form, the equations read 

1^ + V • (pu) = A^VV, (23) 

- 2i7 X u- i^8V^u + ;/2(V^u-f -VV • u), (24) 

3 

^ + u • Ve = -(7 - l)eV • u + K— + A + - pA, (25) 
ot p 
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— = V X (u X B) - V*B + /iV^B, (26) 

V^^-p-l, (27) 

P - (7 - l)pe, (28) 

rd(x,<) = ro, (29) 

{Ti if /9(x,to) > Per 

and < i - to < , (30) 
otherwise, 



where M is the Mach number of the velocity unit, taken equal to unity, J = l/Lj 
is the box size in units of the Jeans length, (j) is the gravitational potential, /2 
is the Galactic disk rotation rate, 1^2 and i^s are respectively the second-order 
and eight-order (hyperviscosity) coefficients, e is the specific internal energy, K 
is the thermal conductivity, and -Tg are respectively the diffuse background 
heating rate and the local stellar heating rate, and /q and Ji are constants. The 
equations are non-dimensionalized to the box size Lq, the velocity unit uq = c, 
the temperature unit Tq = lO^K and the magnetic field unit Bg = 5/xG. 

The cooling rate A is still given by eq. (|]), although the nondimensional 
coefficients are now given in Table |[ together with the non-dimensional heating 
rates, for three different box sizes, at the same temperature unit. 



Table 2. Cooling function parameters for ISM simulations 



Box size (pc) 


C12 


C23 


C34 


C45 


To 


A 


10 


56.2 


0.462 


0.102 


0.474 


4.57 X 10"^ 


20 


150 


845. 


6.95 


1.54 


7.14 


0.688 


38.7 


1000 


5.62 X 10^ 


46.2 


10.2 


47.4 


4.57 


250 



The ISM simulations are started with uncorrelated gaussian random fluctu- 
ations in all variables and amplitudes of order unity, with characteristic scale 
~ 1/8 of the box size. The initial magnetic field has a uniform component of 1.5 
nG on the x-direction, and an rms fluctuation amplitude of 4.5 /iG. The energy 
injection mechanism consists of small-scale 10 pixels across) heat sources 
turned on at sites where the density exceeds a certain threshold (chosen within 
the cold stable branch of the density range), that remain on for At = 6 Myr, 
except in the simulations with box size = 10 pc, for which the time unit is too 
short, and a star would remain on for more than half the duration of the sim- 
ulation. In this case we have shortened At by a factor of 10, and increased the 
energy injection rate by roughly the same factor. These sources are intended to 
mimic the effect of ionization heating from OB stars in HII regions. Supernova- 
like sources are not included because of limitations of the numerical scheme to 
handle very strong shocks. 



26 Vazquez-Semadeni et al. 



6.1 Numerical considerations 

The numerical simulation of thermally unstable turbulent flows presents a signif- 
icant numerical challenge because it requires solving simultaneously the regions 
of sharp gradients occurring in the immediate neighborhood of clouds and the 
quasi- linear development of perturbations in the more distant, relatively quies- 
cent, unstable medium mediating the clouds and the warm, stable, diffuse phase, 
the question being whether this medium fragments, to finally end with a state 
in which the unstable gas virtually disappears from the simulation. The simul- 
taneous solution of both regimes is important, because stellar energy injection 
recycles gas from the dense phase into the warm phase. 

The simulation of regions with strong shocks requires the use of artificial 
viscosities and diffusivities in order to spread out ( "capture" ) shocks over a few 
grid points. Unfortunately, such diffusivities also have the effect of damping the 
growth of perturbations in the relatively smooth regions, because they artificially 



increase the Field length (§2.1), reducing the range of unstable scales, as well as 
their growth rates, in the unstable gas. This problem is also present when finite- 
difference schemes, which produce numerical diffusivity, are used. This was not 
a problem, however, in sections ^ and ^ because of the very large resolutions 
used there, at the expense of using a ID approach. 

The Field length Ap has the unfriendly (for our purposes) property of being 
much larger than the diffusive scale (understood as the molecular mean free 
path) when the latter is very small, since it can be easily shown |Q that the 
wavenumber associated to the Field length satisfies 



(31) 



where fc/f, given by eq. (Q), is the wavenumber associated with the molecular 
mean free path. Thus, fcp only grows as the square root of kx- In a numerical sim- 
ulation, the scale associated with the artificial thermal diffusivity Xk = 2t: /kK 



(cf. §2.1) plays the role of the molecular mean free path, and therefore the growth 
of perturbations with wavelengths Xk < Ap can be artificially suppressed if the 
numerical Ap is much larger than the real Ap in the atomic ISM. In other words, 
much of the resolution in a simulation is "wasted" , in the sense that intermediate- 
wavelength perturbations will be damped even if the diffusive scale is comparable 
to the smallest resolved scale, as we will show below. 

On the other hand, the pseudo-spectral numerical scheme we use in simula- 
tions of the ISM has the advantage that it produces no numerical diffusivities 
at all (the spatial derivatives are calculated exactly in Fourier space, rather 
than through finite-difference approximations). Thus, all artificial diffusivities 
are included explicitely in the equations, and can be controlled through their as- 
sociated coefficients. This allows us to perform simulations of the non-diffusive 
case in the linear and weakly nonlinear regimes (recall the diffusivities are only 
needed to smooth out strong gradients) and of diffusive cases, so that we can 
both test the code against the predictions of the linear analysis, and then mea- 
sure the effect of the diffusivities precisely. We have found that actually the mass 
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difFusion has a much stronger damping effect on the growth of perturbations than 
heat conduction and the viscosity. The value of the thermal diffusivity used in 
the ISM simulations discussed below, K = 6.6 x 10~^, was found to not affect 
the growth rate by more than 10% for the wavenumbers tested below. A similar 
situation holds for the viscous and hyperviscous coefficients. 

We have therefore performed many simulations at low resolution (128 x 128) 
to measure the growth rates of pure sinusoidal isobaric perturbations of various 
wavelengths (from 1 to 1/32 x the box size) as a function of the mass diffu- 
sion coefficient, in order to investigate at what point the perturbation growth is 
suppressed and thus full ISM simulations using those values cannot be trusted 
anymore concerning the instability of regions of sizes comparable to those per- 
turbations. Note that higher resolution is not necessary for this purpose, as all 
that is needed is to resolve the perturbations themselves and their initial growth. 
Higher resolution is only necessary in the fully nonlinear case, to resolve shocks 
while still having a large range of scales between the simulation size and the 
scale of shock-spreading. 

Together with the theoretical growth rate, figure |l] shows the growth rates 
(defined as the inverse of the e-folding time for linear - 2.5% amplitude - per- 
turbations) measured in simulations of a region of 150 pc on a side as a function 
of wavenumber and of the mass diffusion coefficient. These rates can be com- 
pared with the solid line, which gives the solution of the dispersion relation (||) 
as a function of wavenumber. A good agreement is seen between the theoretical 
and numerical rates for the zero-diffusivity case. For non-zero diffusivity, the 
numerical rates are seen to decrease significantly, by roughly a factor of 1.5 at 
/X = 4 X 10^"^ and A = 1/8 of the box size. For the same diffusivity, perturbations 
of size 1/16 of the box size have a zero growth rate (the perturbations remain 
static, without growing or dispersing). 

This damping effect is alleviated somewhat by noting that larger-amplitude 
perturbations have larger growth rates [7^ . We have thus also computed the nu- 
merical growth rates for perturbations of initial amplitude of 25% in simulations 
of box size = 150 pc (not shown), verifying the occurrence of larger growth rates 
in this case. This is advantageous because then weakly nonlinear perturbations 
of scales down to sizes 15 pc should grow in simulations with box sizes 150 
pc at rates comparable to those of linear perturbations without any diffusiv- 
ity. Nevertheless, we have found that zero growth still occurs at the same scale 
(1/16 of the box) as for the linear perturbations at roughly the same value of 
II. Any scales below this are stabilized by the mass diffusion, and thus may be 
wrongly interpreted as stable in a 150-pc ISM simulation using this value of the 
diffusivity. 

Figure ^ shows the numerical growth rates for 2.5% perturbations on simu- 
lations of box sizes 10 pc and 1 kpc on top of the theoretical growth rate curve. 
Remarkably, it is seen that the mass diffusion is much more effective in damping 
the growth of perturbations for a small box size than for a large one. Specifically, 
a diffusivity of /i = 0.001 is enough to completely damp perturbations of size 1/8 
of the box when the box size is 10 pc, while a value fi = 0.008 still allows growth 



28 Vazquez-Semadeni et al. 



of perturbations of size 1/16 of the box for a box size of 1 kpc. This is actually 
easy to understand because a larger box size is represented, in non-dimensional 
units, by larger thermal coefficients (see Table |^), while the diffusive coefficients 
are independent of the physical box size, if the diffusive scale is kept at a given 
fraction of the numerical box (in pixels). This is similar to the effect on the 
Field length of the thermal conductivity (eq. ||3^]). Somewhat surprisingly, we 
conclude that, at a fixed value of the mass diffusion coefficient, a simulation with 
a larger physical box size is more accurate than one with a small physical size. 

A final comment is that, from the results of this section, it appears that 
the numerical scheme best suited for solving the problem at hand (i.e., captur- 
ing strong gradients occurring in dense clouds while not disturbing the growth 
of perturbations in the mildly turbulent diffuse medium surrounding the dense 
clouds), is a pseudo-spectral scheme with a position-dependent value of the dif- 
fusivities, as in eq. (|l6|), so that strictly zero-diffusivity can be used in mild 
regions, while still smoothing out the strong gradients in the clouds can be 
achieved. Wc will attempt such an approach in forthcoming papers. However, 
our presently available tools for modeling the full ISM problem only include 
constant-coefficient diffusivities, and we will discuss results using them in the 
next subsection. Note that not even adaptive-mesh refinement schemes may be 
better suited for this problem, because the small-amplitude perturbations need 
to be followed in the diffuse medium in order to determine whether they grow 
spontaneously, and this would require high refinement levels in relatively large 
volumes. 



0.5 F 




0.001 0.010 0.100 1.000 10.000 100.000 1000.000 
1/l[pc] 

Fig. 10. Numerical growth rates for simulations of the linear growth of 2.5% isobaric 
density perturbations in boxes of sizes 10 and 1000 pc and various values of the mass 
diffusion coefficient, superposed on the growth rate curve. Crosses: 1000 pc, fi — 0. 
Asterisks: Box size=1000 pc, ^ = 0.008. Triangles: Box size=10 pc, ^ = 0.0. Diamonds: 
Box size=10 pc, fi = 0.001. Squares: Box size=10 pc, fj, = 0.0005. 
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6.2 Results 



In the light of the results of i: 6.1, we can now proceed to present some numerical 
simulations of the ISM and assess their reliability regarding the mass fraction in 
the unstable range. 

In Paper IV, we presented a 2D simulation of the ISM with a 1-kpc box at a 
resolution of 512^ grid points with a mass diffusion coefficient /i — 0.0075. From 
fig. |l^ it can be seen that perturbations of sizes down to 1/16 of the box can 
grow at this value of [i with a physical box size of 1 kpc, while perturbations of 
size 1/32 of the box are damped. Thus, although this box size is the one that 
provides the largest range of unstable scales, it nevertheless does not reach down 
to the fastest-growing scales. Thus, we have performed two more simulations, to 
cover the entire range of scales of interest: the first one with a box size of 150 
pc, and a very-high resolution of 1536^ grid points, in order to be able to use a 
mass diffusion coefficient [l — 0.003 which, according to the data in fig. |], allows 
growth again of perturbations of size 1/16 of the box, thus barely reaching the 
scales of fastest growth ( ;S 10 pc). The second simulation uses a box of 10 pc, and 
/U = 5 X 10~^ at a resolution of 512^, thus allowing the growth of perturbations 
down to scales 1/8 of the box size. We have opted for performing 2D simulations 
in order to maintain relatively high resolutions, while still being able to capture 
the complex vector interactions of the system. 

In fig. |ll] we show the density and temperature histograms of these three 
runs after a stationary regime has been attained. It can be seen that the density 
histograms do show signatures of TI, such as slope changes and slight peaks 
at the stable densities, but nevertheless, a sizeable fraction of the gas is in the 
unstable regime. A similar result was reported by Kritsuk & Norman p6[ . 

Concerning the temperature histograms, we see that, in fact, as the physical 
box size decreases, the temperature histogram has less pronounced spikes at the 
temperatures corresponding to the stable phases, suggesting that phase segre- 
gation is less pronounced as well. In the case of the 10-pc run, this can be an 
artifact of the reduced unstable range due to the mass diffusion. However, this 
is not so for the 1-kpc and 150-pc runs, as both have roughly the same range of 
unstable scales, and in fact with larger growth rates in the case of the latter, be- 
cause the corresponding physical scales are smaller. Thus, these two runs suggest 
that in fact the result is real, due to the decreasing value of ry as smaller scales 
are considered, because in this case the perturbations are adiabatic-like and are 
stable to first order. Taking the simulations at face value, the mass fractions in 
the "unstable" temperature range are 27% for the 1-kpc run, 58% for the 150-pc 
run, and 92% for the 10 pc run. Again, we see the trend of a larger fraction of 
unstable gas as smaller scales are considered. 

On the other hand, the effect of strong compressions continues to promote the 
instability. Figures |lj and |l^ respectively show the density and pressure fields for 
the 150-pc run at a typical time. There it can be appreciated that the density 
maxima (and the filaments connecting them) correspond to pressure minima^ 
indicating that these regions are unstable, even when they are much smaller 
than the smallest linearly unstable scale allowed by the mass diffusion. This 
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Fig. 11. Top: Density histograms for the ISM simulations with box sizes 10 pc [solid 
line), 150 pc [dotted line) and 1000 pc [dashed line). A clear trend toward higher phase 
segregation (more strongly bimodal shape) is seen at large physical box sizes because 
the wave mode is unstable at large scales only. Bottom: Temperature histograms, with 
the same labeling. 



is because the strong turbulent compressions locally increase the density and 
the cooling rate, decreasing rj. Thus, these regions belong to the long-wavelength 
(small-rj) regime in spite of having small physical sizes, and the pressure behaves 
closer to the thermal-equilibrium curve. These small-scale, small-77 regions are 
analogous to the compression-induced instability of Hennebelle & Perault 
except that in this case they have been pushed from the large- to the small-ry 
regime by the compression, rather than from the stable to the unstable regime. 
They also correspond to the transition from isobaric to isochoric cooling as the 
density increases described by Burkert & Lin ||] and Kritsuk & Norman |46 



It is worth comparing here the present results to those of Hennebelle & 
Perault in somewhat greater detail, as those authors did not find significant 
amounts of unstable gas in their simulations after the condensation process 
ended. This is because their setup did not consider a globally turbulent medium, 
but only the effect of a single compressive wave of intermediate-strength (Mach 
number ~ 2) on the already-segregated warm stable phase. Instead, here we are 
considering the case of a medium with mean density n ~ 1 cm~^, which is close 
to the ISM average density and lies in the unstable range, so that, even if the 
two phases are segregated, the mean density remains at that value. Moreover, 
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rather than the effect of a single large-scale compressive mode, we are consid- 
ering a globally turbulent medium in which the source of energy is stellar-like, 
having the effect of recycling matter from the dense phase into the diffuse one. 
As shown in the simulations, a situation similar to that of Hennebelle & Perault 
applies at the sites where the compressions are strong (near stellar sources), ex- 
cept that in their case the triggering is in absolutely linearly stable gas, while 
in our case the triggering is in gas stable only to adiabatic perturbations. On 
the other hand, at more remote sites, where the turbulence is weaker the 
fluctuations remain stable in general. 




Fig. 12. Density field of the 150-pc simulation of the ISM at the same time for which 
the density and temperature histograms are shown in fig. ^} The resolution is ISSG'^ 
grid points. 

It should be emphasized that, in the context of full ISM simulations, the 
regions with a reversed pressure gradient are seen to occupy a very small fraction 
of the volume, with the majority of the space being occupied by a moderately 
turbulent medium in a nearly isobaric regime, but with significant fluctuations 
around it. Particularly noteworthy is the existence of weak expanding shock 
waves which propagate away from the filamentary clouds, behind which both the 
density and the pressure are slightly increased with respect to the intercloud 
medium, similarly to the case of the wave observed in run DEN75 in §|. It 
should be noted, however, that the origin of the outgoing shock waves is not 
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Fig. 13. Thermal pressure field of the 150-pc ISM simulation. The highest pressures 
correspond to regions of "star formation" clustering, in which the heating from many 
stars combines additively. Note the high-pressure regions around the clouds, left by the 
weak shocks propagating into the intercloud medium. 
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Fig. 14. Pressure vs. density for a vertical stripe of width 200 pixels (~ 20 pc) in the 
150-pc simulation, passing through the large complex slightly up and to the left of the 
simulation center (cf. fig. [L2[). 



Thermal Instability and Magnetic Pressure in the Turbulent ISM 



33 



exactly the same in the ISM simulations as it was in run DEN75. In the latter, 
the shock originates when the condensation reaches its maximum (overshooting) 
density and terminates the converging motions within it. In the ISM simulations, 
it originates when star formation suddenly heats and repressurizes the clouds. 
Nevertheless, the effect on the gas surrounding the clouds is similar. The volume 
between the wave front and the central clouds contains gas with a "regular" 
pressure gradient (S^) and temperatures in the "unstable" range. However, as 
in the case of run DEN75, this gas is not truly unstable, as it is out of thermal 
equilibrium, albeit close to pressure equilibrium. The existence of gas with a 
regular pressure gradient is seen in fig. which shows the thermal pressure vs. 
the density for a vertical stripe of width 200 pixels 20 pc) and length equal to 
the box size in the 150-pc simulation including most of the large complex slightly 
up and to the left of the simulation center. It can be seen that the pressure is far 
from having a unique value, and, in particular, has a scatter of about one order 
of magnitude in the "unstable" range, with the upper envelope of the points 
having a regular dependence of pressure on density. This effect is even more 
pronounced in the simulations of Mac Low et al. [ pT| , which include supernova 
energy input (see also the chapter by Mac Low, this volume). 

In order to complete the description of the dynamics of the low-thermal- 
pressure (Iow-Pt) regions, it is important to also analyze the magnetic pressure 
within them. This is shown in fig. Interestingly, some of the \ow-Pt filaments 
are also seen as low-magnetic-pressure filaments, while some others are not. The 
cases with reduced magnetic pressure can be interpreted as sites of field reversals, 
which arise when the compressions amplify the (small-scale) field fluctuations 
perpendicular to the direction of compression, as observed by Passot, Vazquez- 
Semadeni & Pouquet |66j. As the field reverses, it passes through zero. However, 
very few examples of a filament seen also as a high-magnetic-pressure structure 
are seen, although the magnetic pressure maxima are indeed seen to occur in 
the star-forming regions. We thus conclude that in general the collapse of the 
Iow-Pt regions is ensured. This is confirmed by the fact that those filaments 
systematically become sites of new star formation at later times as the density 
threshold is reached. 

Also of interest is that the enhanced-density, enhanced- regions surround- 
ing the clouds in general also have an enhanced magnetic pressure, so that 
these regions are in general slightly over-pressured with respect to the global 
intercloud medium, and expand somewhat. The intercloud medium is com- 
pletely permeated by these traveling fronts originating from the clouds. This 
is most clearly seen in animations, which can be retrieved from our web site 
jittp : //www ■ astrosmo . unam . mx/^e . vazquez/turbulence_HP . Finally, the mag- 
netic field is rather uniform in the intercloud medium. More discussion concern- 
ing the magnetic pressure in ISM simulations can be found in the chapter by 
Mac Low. 
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Fig. 15. Magnetic pressure field of the 150-pc ISM simulation. Note that some filaments 
of low thermal pressure (fig. also have low magnetic pressure, while some others 
have undisturbed magnetic pressure. 

6.3 Discussion and caveats 

The results from the previous section suggest a dynamic ISM in which several 
physical processes are occurring simultaneously besides the pure condensation of 
thermally-unstable regions, and in which, in fact, the very nonlinear development 
of the latter leads to complex dynamics. 

These results are not free of caveats, however. Three main limitations of the 
simulations prevent the results from being definitive. First, the need to use artifi- 
cial diffusivities, especially in the continuity equation, to broaden discontinuities 
out to a few grid points, enormously shortens the range of unstable scales, typ- 
ically stabilizing scales from the resolution limit up to 1/32 or 1/16 of the box. 
This implies that the fraction of unstable gas could possibly be overestimated 
because scales that could fragment in the real ISM do not in the simulations. 
However, "unstable" structures larger than those damped by the mass diffusion 
are seen to exist in the simulations as well, suggesting that the effect is real, if 
perhaps not as strong. 

Here it is important to emphasize that at large scales the medium is in- 
deed unstable, but this is manifested not in the spontaneous condensation of 
the structures, but in their null resistance to compression, since typically the 
turbulence velocity dispersion is supersonic, and the structures are compressed 
by the turbulence before they can spontaneously condense. At small scales, on 
the other hand, as indicated by the runs in a moderate amount of turbulence 
prevents the spontaneous condensation because the fluctuations are closer to 
being adiabatic than isobaric. 

A second caveat is that the threshold density for star formation (SF) used in 
the simulations is rather low (ricr = 15 cm^"^ for the 1-kpc and the 150-pc runs. 
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and 7icr ~ 25 cm^'^ for the 10-pc run), again in order to avoid extreme gradients 
that would form upon the onset of stellar heating in very smaU, dense clouds. 
This implies that the clouds are not allowed to reach the pressure-equilibrium 
density in the dense phase. The SF scheme thus bypasses the process of a cloud 
(or cloud complex) becoming self-gravitating to start forming stars. This presents 
the risk that perhaps not enough accumulation of mass is allowed in the clouds 
before they engage in SF, forcing a greater fraction of the gas to remain in transit 
in the unstable range towards the clouds. We feel, however, that the essence of 
the process is the fact that matter is recycled continually among the phases 
and, given the time scales of the problem, that the existence of sizeable amounts 
of gas always traversing the unstable range is inevitable to some extent. More 
accurate estimates of this fraction should be possible upon the application of 
modifications to the numerical method, as decribed above. 

The third caveat is the omission of supernova- (SN-)like energy imput, which 
should sweep a larger fraction of the volume and induce the formation of a 
hot phase, which we have neglected. Their inclusion would cause that a larger 
fraction of the volume would be swept up into shells, forcing the conversion of 
diffuse gas into the dense phase. Nevertheless, the turbulent mixing ought to 
be even stronger in this case, and relatively quiescent regions must still exist, 
as observations suggest. These would undergo slow condensation in the manner 
outlined here. 

In any case, the results presented here should only be considered as sugges- 
tive of the existence of sizeable amounts of unstable gas in the ISM, and more 
specifically designed numerical methods should be applied in order to obtain 
more definitive confirmation. 

7 Summary and conclusions 

In this paper we have reviewed results from various studies aiming at under- 
standing the role of TI in the turbulent atomic ISM, and the behavior of the 
magnetic pressure in the fully turbulent case. The motivation has been twofold. 
On the one hand, the classic multi-phase models of the ISM have neglected the 
implications of the ISM being turbulent, and it is thus important to assess the 
consequences of advection on the thermal and spatial structure of this medium. 
On the other hand, observations have often suggested the presence of gas with 
temperatures in the thermally unstable range, in apparent contradiction with 
the multi-phase models. 

We first reviewed the classic instability analysis of Field emphasizing 
the different behavior of long- and short-wavelength perturbations (for which 
the ratio rj of the cooling [tc] to the sound crossing [ts] time is respectively small 
and large), and of entropy and isentropic perturbations (which trigger the con- 
densation and the wave modes, respectively). We pointed out that, while much 
study has been devoted to isobaric entropy perturbations, real- world fluctuations 
in the ISM are produced through velocity fluctuations which, in the small-scale 
limit, belong to the isentropic kind, and are therefore stable to first order at 
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small scales. We also briefly reviewed the magnetic case, in which the presence 
of a uniform magnetic field can stabilize perturbations with wavcmimbers per- 
pendicular to it. 

We then reviewed results on the nonlinear stages of evolution of isobaric 
entropy perturbations, focusing on those that have quantified the magnitude of 
the speeds developed and the times required for completing the condensation 
process as a function of the parameter ?/. For presently-accepted values of the 
heating and cooling rates [^, large-scale initial perturbations ( ^ 15 pc, rj ^ 0.2) 
develop supersonic speeds, require times <L 10 Myr to complete the condensation 
process, and end up with densities and pressures above the thermal equilibrium 
value due to the ram pressure of the still infalling gas. Those times are long 
compared with typical times between successive external shock passages, and star 
formation time scales. We thus concluded that clouds formed from perturbations 
of such sizes (although the resulting cloud has a size ^ 1 pc) are unlikely to exist 
in thermal pressure equilibrium with their surroundings. Initial perturbations 
of sizes ^ 3 pc, on the other hand, require times ~ 4 Myr to complete their 
evolution and do not generate supersonic speeds, thus reaching a more quiescent 
final state, and adhering better to the paradigm of thermal-pressure bounded 
clouds at the end of their evolution, although, by the time the cloud has formed, 
accretion is still occurring, so that the clouds are bounded by weak accretion 
fronts, rather than contact discontinuities. Furthermore, the gas still accreting 
is necessarily in the unstable temperature range, although it is not in thermal 
equilibrium; instead, it has a "regular" pressure behavior (P increases with p), 
and thus it is not prone to further fragmentation. 

We then described the evolution of perturbations induced by turbulent ran- 
dom forcing. In this case the crossing time entering 77 should be taken as the 
minimum of the sound and the turbulent (tu) crossing times. Thus, small-scale 

0.3 pc) velocity fluctuations are quasi-isobaric at very small amplitudes, be- 
cause in this case > Tq > Tg so that the flow can cool in response to the velocity 
perturbation. As the perturbation amplitude is increased, so that Tc > Tg > Tu, 
the situation changes because now the density is driven by the turbulent veloc- 
ity rather than by sound waves, and the perturbations become quasi-adiabatic 
in character, becoming stable. We empirically found this to occur roughly when 
the rms Mach number ^ 0.3. Finally, however, if the perturbation amplitude 
becomes very large, then the density increment induced by it becomes nonlinear 
and accelerates the cooling rate, effectively causing 77 < 1. In this case, velocity 
fluctuations trigger the condensation mode, which is again unstable, and cause 
condensations. 

Thus, we reached the important conclusion that small-scale fluctuations be- 
have very differently when they are entropy perturbations (caused, for example, 
by local variations in the heating or cooling rates) and when they are adiabatic 
(caused by velocity fluctuations), being unstable (and with the fastest growth 
rates) in the former case, but linearly stable in the latter. 

We then considered the magnetic field as an additional source of pressure in 
the ISM, confirming earlier results that at low and intermediate densities the 
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magnetic pressure is strongly decorrelated from density in fully turbulent cases 
(large field fluctuations), and proposed an interpretation of this phenomenon in 
terms of the scaling of with density for the slow and fast modes of simple 
nonlinear MHD waves. The decorrelation between magnetic pressure and den- 
sity has several implications, among which is that the magnetic field probably is 
ineffective in supplementing thermal pressure in highly turbulent, thermally un- 
stable conditions, and that it is probably inadequate to model magnetic pressure 
by means of an equivalent polytropic behavior in the fully turbulent case. 

Finally, we discussed results from simulations of the ISM in more than one 
dimension at large and intermediate scales and at various resolutions. To this 
end, we first performed a detailed study of the competition between numerical 
diffusivities and the growth of TI, finding that even when the diffusivities (es- 
pecially the mass diffusion, which is necessary numerically) are confined to the 
smallest scales on the numerical grid, they can push the smallest unstable scale 
(the "Field" length Af) to relatively large scales in the simulations, especially 

1/2 

for small physical simulation sizes, because Ap scales as A^ , where Xk is the 
diffusive scale. 

With this information, we discussed the fact that many ISM simulations 
suggest that the basic structure does not depend sensitively on whether TI is 
present, as long as there are turbulent motions driven by stellar- like sources (that 
imply recycling of gas from the cold to the warm phase), and that significant 
fractions of the gas mass (15-50%) appear to be in the unstable regime. This 
appears to be a consequence of the fact that the diffuse medium is in a moder- 
ately turbulent state, so that a) the fluctuations there have a regular pressure 
gradient and b) the magnetic field is not strongly turbulent, and therefore may 
cause additional stability. Of course, when the relatively quiescent intercloud 
medium is hit by a strong shock from, say, a supernova remnant, then TI can be 
rapidly induced, as in the studies by Hennebelle & Perault |^2|,^ and Koyama 
& Inutsuka Q. 

A final remark of interest is that it may be possible to determine obser- 
vationally whether the gas seen at unstable temperatures corresponds to the 
out-of-thermal-equilibrium gas observed in the simulations by either a) simulta- 
neously determing two of its thermodynamic variables, or b) comparing directly 
observed cooling rates (e.g., fine structure hues) with theoretical estimates of the 
heating rate (e.g., photoelectric heating) in specific regions (C. Heiles, private 
communication). If this is confirmed, then it would provide strong evidence that 
turbulent motions populate all regions of the thermodynamic variable space, 
preventing a sharp segregation of the atomic ISM into the stable phases of TI. 

We have greatly benefitted from exchanges with C. Heiles, P. Hennebelle, 
H. Koyama, J. Scalo and E. Zweibel. The report from an anonymous referee 
prompted much improvement of the paper and led us to the study of numerical 
damping of the growth rates. This work has received partial financial support 
from CONACYT grant 27752-E, from the French national program PCMI, and 
from the conference organizers to E.V.-S. We have made extensive use of NASA's 
Astrophysics Data System Abstract Service. 
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